Expanding the limits of nuclear stability at finite temperature

Properties of nuclei in hot stellar environments such as supernovae or neutron star mergers are largely unexplored. Since it is poorly understood how many protons and neutrons can be bound together in hot nuclei, we investigate the limits of nuclear existence (drip lines) at finite temperature. Here, we present mapping of nuclear drip lines at temperatures up to around 20 billion kelvins using the relativistic energy density functional theory (REDF), including treatment of thermal scattering of nucleons in the continuum. With extensive computational effort, the drip lines are determined using several REDFs with different underlying interactions, demonstrating considerable alterations of the neutron drip line with temperature increase, especially near the magic numbers. At temperatures T ≲ 12 billion kelvins, the interplay between the properties of nuclear effective interaction, pairing, and temperature effects determines the nuclear binding. At higher temperatures, we find a surprizing result that the total number of bound nuclei increases with temperature due to thermal shell quenching. Our findings provide insight into nuclear landscape for hot nuclei, revealing that the nuclear drip lines should be viewed as limits that change dynamically with temperature.


Supplementary Note 1. THEORETICAL FORMALISM
We consider the nucleus at finite temperature within the relativistic energy density functional (EDF) formalism. The most common approach at finite temperatures is to solve the finite-temperature relativistic mean-field equations (FT-RMF). However, with non-vanishing temperatures, a proper treatment of the continuum states is mandatory to achieve convergence of our results. Therefore, we implement the method first introduced by Bonche, Levit and Vautherin in Refs. [1,2], further denoted as the BLV prescription. Finally, we generalize the BLV prescription to treat superfluid nuclei within the finite-temperature relativistic Hartree-Bogoliubov (FT-RHB) theory, including also the deformation and pairing in the calculations.

A. The FT-RMF theory
The FT-RMF equations are derived by assuming the nucleus is an isolated system that cannot exchange heat or particles with the environment [3]. Such a system is characterized by its temperature T and chemical potential λ. In order to derive the FT-RMF equations, we start by defining the grand-canonical potential Ω where the internal energy is E = Tr DĤ , entropy is S = −k B Tr D lnD and particle number is N = Tr DN ,D being the single-particle density operator and k B is the Boltzmann constant. The chemical potential is defined both for protons (λ p ) and neutrons (λ n ). The many-body Hamiltonian operator assumes the usual form [4] where t ij is the kinetic energy andv ijkl = v ijkl −v ijlk anti-symmetrized two-body interaction, while c † k , c k are the single-particle creation and annihilation operators, respectively. The gist of the FT-RMF theory is to assume an ensemble of independent-particles; thus, the density operator is approximated by the FT-RMF density operator [3] where the Fermi-Dirac factor is defined as f i = [1 + exp [β(ε i − λ q )]] −1 , β = 1/(k B T ), while ε i is the single-particle energy. The number operator of a particular state is defined aŝ which is also known as the generalized Wick contraction [5]. By evaluating the respective traces within the FT-RMF approximation, we can rewrite the grand-canonical potential as with the usual notation Γ ij = klv ikjl ρ lk . We perform variations over the density ρ and Fermi-Dirac factors f to get [3] δΩ = Tr [(t − λ q )δρ] + Tr [Γδρ] with the mean-field Hamiltonian defined as h = t + Γ − λ q . Within the RMF formalism above equation takes a form of the single-particle Dirac equation [6,7] [−iα∇ + βM * (r) + V (r)] ψ i (r) = ε i ψ i (r), where the effective mass is defined as M * (r) = m + S(r), S(r) is the scalar potential, V (r) the vector potential and m the bare nucleon mass. The single-particle wavefunction where s denotes the spin coordinate and Φ n (r, s) is the harmonic oscillator wavefunction. The expansion in upper(lower) components is cut-off at some shell number n max (ñ max ).
where B (j) n (r, s) is a generalized interpolation function of j-th order.
For now, we keep the equations general and later we will assume a specific symmetry of a nucleus: either spherical or axially-deformed. Details of the vector V and scalar S fields depend on the selected functional. In this work, it is either the meson-exchange (DD-ME2) or the point-coupling (DD-PC family).

B. Continuum subtraction -the BLV prescription
As the temperature increases, a non-vanishing number of single-particle states gets scattered to the particle continuum. In our work, we employ the prescription developed in Refs. [1,2], which we denote as the BLV prescription in the following. Within the BLV prescription, we can define two density matrices: (i) ρ density matrix, associated with the system consisting of a nucleus and a surrounding vapor (Nucl+Vap), and (ii)ρ which corresponds to the vapor solution (Vap) only. Minimization is performed with the subtracted grand-canonical potential defined as ∆Ω = Ω(ρ) − Ω(ρ). (S11) By inserting Eq. (S5), we obtain: where we define the subtracted density asρ = ρ −ρ. The subtracted entropy is given bȳ (S13) The summation is performed over the discrete single-particle states (denoted by k) and integration is over the continuum states with level density g(ε). However, in practical calculations, continuum is often discretized, therefore we replace the integration with summation over the continuum states. Using the short-hand notation of Ref. [4], we rewrite the second term in Eq. (S12) as 1 2 Tr Γρ −Γρ = 1 2 Tr 1 Tr 1 (ρvρ) + Tr 1 Tr 1 (ρvρ). (S14) The chemical potential λ q is fixed in order to reproduce the total particle number N q First, to get the FT-RMF equation for the Nucl+Vap system, we vary Eq. (S12) with respect to the Nucl+Vap density ρ and the FT-RMF equation for the Vap system is obtained by variation with respect toρ The coupling between the densities of these two systems is given through Eq. (S15). However, in Refs. [1,2] it is noted that the Coulomb term has to be treated separately in order to avoid divergences, stemming from the long-range charged vapor contribution. The BLV prescription imposes the following form of the Coulomb interaction term (S18) Components of the nuclear interaction remain unchanged. For demonstration purposes, in the following, we assume that the only interaction term at the mean-field level is the direct Coulomb term (and denote it by v c with the corresponding proton density ρ p ). Varying Eq. (S12) with ρ andρ gives us the FT-RMF equations for the Nucl+Vap and Vap system, (S19) Therefore, within the BLV subtraction procedure, the self-consistent FT-RMF problem con- The pairing interaction can be included in a straightforward manner by introducing the Bogoliubov transformation of single-particle operators [4] where β † k , β k are the quasi-particle (q.p.) creation and annihilation operators, respectively. We can now introduce the pairing tensor κ, which together with the single-particle density in the single-particle basis {c † k , c k } assumes the form [3] ρ By introducing the pairing correlations, Eq. (S12) can be rewritten as where κ,κ is the pairing tensor for the Nucl+Vap and Vap systems, respectively, and ∆ ij = where (U, V ) is the q.p. wave function and E the corresponding q.p. energy. On the other hand, by varying the Eq. (S22) withρ andκ we get the FT-RHB equation for the vapor where (Ũ ,Ṽ ) is the vapor system wave function andẼ is the eigenvalue. For the pairing part of the nuclear Hamiltonian we employ the separable interaction, first introduced in Ref. [8], of the form where R = 1/2(r 1 + r 2 ) and r = r 1 − r 2 denote the center-of-mass and relative coordinate, respectively, while P (r) has the form P (r) = 1 The parameters G and a used together with the DD-ME2 and DD-PC1 interactions are taken from Ref. [8], while the parameters used with the DD-PCX interaction are from [9].
For the DD-PCJ family of functionals those values can be found in Ref. [10].
D. Comparison of the relativistic and non-relativistic models without the BLV subtraction We first compare the relativistic and non-relativistic model calculations using the FT-RHB and finite temperature Hartree-Fock-Bogoliubov (FT-HFB) methods without the BLV procedure. For comparison, the results for the Skyrme-type SkM* interaction are taken from Ref. [11]. Since the FT-HFB results shown in Fig. 8 in Ref. [11] assume spherical symmetry, the FT-RHB calculations are also performed by assuming spherical symmetry to ensure the fairness of the comparison. In Figure 1, the two-neutron separation energies (S 2n ) and the corresponding chemical potentials (λ n ) are displayed for the selected isotopic chains with increasing neutron number and temperature using SkM* (a)-(b) and DD-PC1 nuclei with λ n ≤ 0, which is consistent with the approach in Ref. [11]. It is clear that both relativistic and non-relativistic models predict the same behaviour with increasing neuron number and temperature. At zero temperature, the S 2n and λ n increase with increasing neutron number, as expected. The sharp changes in the S 2n and λ n values can be seen around the nuclear magic numbers. By increasing temperature, the shell effects start to disappear, and the changes in the S 2n and λ n values become smoother. At high temperatures, i.e., at T = 2.0 MeV, the number of nuclei with negative chemical potential or bound nuclei increases. On the other hand, the predictions using the S 2n values are not consistent with the behaviour of chemical potentials, especially for neutron-rich nuclei near the drip lines.
More nuclei become bound when the chemical potentials are considered, but these nuclei also have positive S 2n values and are unstable against two-neutron emission. It is known that continuum effects play an important role in the determination of the position of the drip lines already at zero temperature. At high temperatures, i.e., T > 1 MeV, the impact of the states in the continuum increases due to the nucleons scattered throughout this region.
Therefore, a special caution must be taken in the calculations for weakly-bound nuclei at high temperature. In Supplementary Note 2 C, we show that the inconsistencies between the location of the neutron drip lines using the chemical potential and two-neutron separation energies can be resolved by the inclusion of the BLV prescription in the model calculations, and using the Free energies to calculate the two-neutron separation energies.

Supplementary Note 2. NUMERICAL METHODS AND TESTS
Most of the calculations in this work are performed by assuming the axially-deformed reflection symmetric nuclei. However, since the axially-deformed calculations are very timeconsuming (especially with the explicit treatment of continuum), we initially check our implementations by assuming the spherically-symmetric mean-field at finite temperature.
Two FT-RMF solvers that we employ are: • FT-RMFBSPL -finite-element method (FEM) algorithm which solves the FT-RMF equations by discretizing the wave functions on a B-spline mesh in the coordinatespace. Our code is based on the method presented in Ref. [12]. The basis parameters are: number of finite-elements N fe , order of the B-splines N ord and the box-size R box . • FT-DIRHBz -solves the axially-deformed FT-RHB equations by assuming reflection symmetry. It employs expansion of the wave functions and fields in the axiallydeformed HO basis. The code is based on Ref. [13].
All calculations in this section employ the relativistic DD-ME2 interaction.

A. Convergence with the increasing box-size
For weakly-bound nuclei, even low temperatures can cause results not to converge with the box-size R box due to the scattering of single-particle states to the continuum. Therefore, we have implemented the BLV subtraction procedure in our calculations. To test the convergence of main observables with increasing R box we select a nucleus in the vicinity of the drip-line, 202 Sm, at T = 1.0 MeV. Our finite-temperature calculations in axial-geometry indicate that 202 Sm is in a spherical state at T = 1 MeV, with vanishing pairing correlations.
We perform calculations with the FT-RMFBSPL code using the DD-ME2 interaction. Results are shown in Fig. 2 for the dependence of entropy S (left panel) and the RMS neutron radius ⟨r 2 n ⟩ (right panel) on box-size R box . This is a drastic example that demonstrates the importance of proper treatment of continuum in drip-line nuclei at finite temperature.
In the case where we include the continuum subtraction within the BLV prescription (blue circles), results are independent of the box-size, starting from R box ≈ 20 fm. If the continuum is not subtracted (red squares), we observe an almost linear increase in entropy and an exponential increase in the RMS neutron radius with increasing box-size.
These results can be readily understood when observing the vector density ρ v . In Fig.   3 we display the radial dependence of the total vector density for R box = 30 and 40 fm.
Imposed box-boundary conditions necessitate that the density has to vanish at R box . Due to the contribution of continuum states, without BLV subtraction, one obtains a long tail in density (red dashed line), constant up to R box . As the box-size increases so does the density tail. The single-particle states in the continuum, responsible for the long density tail, behave as a vapor that surrounds the nucleus. Such a vapor is precisely the reason for diverging results shown in Fig. 2. Using the BLV procedure allows one to isolate the vapor (Vap) contribution (blue dotted line) from a combined solution of nucleus and vapor (Nucl+Vap). By subtracting the Vap density from the Nucl+Vap density, we obtain the subtracted density (black solid line) that is independent of the box-size.

B. Comparison between the FT-RMFHO and FT-RMFBSPL solvers
It is well known that the HO basis, due to its Gaussian tail, cannot reproduce the proper asymptotic behavior of the nuclear density [14]. Such a drawback does not pertain to serious issues when studying the nuclei in the vicinity of stability valley, however, it can lead to wrong results for phenomena in weakly-bound nuclei near the drip lines. Such an issue could be resolved by: (i) increasing the number of oscillator shells N osc , or (ii) a local-scale transformation [15]. However, the computational time significantly increases with N osc , while the local-scale approximation has only been applied at zero-temperature for the relativistic EDFs [16]. Even at zero temperature, the transformed HO (THO) basis has significant difficulties for the relativistic functionals, due to contribution of anti-particle states. Later, we will explicitly demonstrate the issue with anti-particles in determining the optimal basis parameters for HO expansion. Therefore, in this subsection we compare the results for important observables at finite temperature obtained with the FT-RMFHO solver that employs the HO expansion to the coordinate-space FT-RMFBSPL solver.
First, we have to select an optimal number of shells N osc and oscillator length b 0 for our comparison. In the case of the HO based solver, the box-size is approximately given by . In the finite-temperature calculations, if the number of shells N osc is fixed, b 0 should be determined as to minimize the subtracted free energyF . Note that the subtracted free-energy is defined as where F N ucl+V ap (F V ap ) is the free energy of the Nucl+Vap(Vap) system. In Figure.  C.

Test of chemical potential
At finite temperatures, the binding energy E(Z, N ) should be replaced by the free energy F (Z, N ), therefore, the two-neutron and two-proton drip lines are defined by respectively. In the BLV method, the subtracted free-energyF is employed. This quantity does not take into account the vapor contribution and can be expressed asF (Z, N ) = E(Z, N ) − TS(Z, N ), whereĒ andS represent the subtracted total binding energy and entropy, respectively. Figure 6 shows the impact of BLV subtraction on the total energy  Isotopic dependence of two-neutron separation energy S 2n (a) and neutron chemical potential λ n (b) for the selected even-even nuclei using using the FT-RHB+BLV method with DD-PC1 interaction.

Supplementary Note 3. CALCULATION OF NEUTRON EMISSION LIFETIMES
Before presenting the results of neutron emission lifetimes, it is important to benchmark our model against existing theoretical calculations. This is especially relevant to this work because we use the harmonic oscillator basis, which does not reproduce the exponential tail decrease of the density at large distances from the nucleus. To achieve this aim, we calculate the neutron emission lifetimes at finite temperatures using two numerical approaches: (i) discretizing Eqs. (S23) and (S24) in the basis of the harmonic oscillator (FT-RMFHO) and (ii) discretizing in the B-spline basis, corresponding to the solution in the coordinate space (FT-RMFBSPL). Calculations in the coordinate-space basis are known for reproducing the tails of particle density, which is of special significance for the study of weakly-bound systems and neutron halos [14,[17][18][19]. Unlike being in the ground state at zero temperature, nuclei at finite temperature are found in highly-excited meta-stable states that can decay either by particle emission, provided that the excitation energy is above the particle-decay threshold, or by gamma emission. In particular, as the temperature increases, more and more nuclei gain a finite width for neutron emission. The neutron emission width Γ n can be obtained from the nucleosynthesis formula [2] Γ n h = n gas ⟨σv⟩, where σ is the neutron capture cross-section, ⟨v⟩ is the average velocity of particles in the external nucleon gas, and n gas is the neutron vapor density calculated as number of neutrons in the vapor divided by the discretization volume. The neutron cross-section can be approximated by the geometric cross section σ = πR 2 [2], where the RMS-radius R is obtained from the FT-RHB calculations. Finally, the neutron emission lifetime can be calculated as τ n =h/Γ n . The statistical velocity is calculated from the finite-temperature canonical single-particle neutron energies ε n assuming the Fermi-Dirac distribution of neutrons f (ε n ), therefore [20] ⟨v⟩ where m n is the neutron mass. We obtain the canonical single-particle neutron vapor states by diagonalizing the neutron vapor densityρ n and transforming the corresponding quasi- particle energies in this basis. Since we lose the notion of independent quasi-particles at finite-temperature, this procedure is an approximation.
To check if the HO basis provides reasonable results in comparison to coordinate-spacebased methods, we choose large temperatures of T = 2 and 3 MeV and assume that all nuclei are spherical and in a normal state (no pairing correlations), which allows for detailed calculations. As we will demonstrate in Supplementary Note 4 such an assumption is reasonable. Reducing the basis size due to spherical degeneracy is especially significant for the coordinate-space code, where the calculation time becomes prohibiting for systematic deformed calculations. High temperatures are chosen to exaggerate the influence of the particle continuum. Calculations are performed for the samarium (Z = 62) isotopic chain.
The FT-RMFHO code is characterized by the number of oscillator shells N osc in which the single-particle basis is expanded, while the wave functions of the coordinate-space code FT-RMFBSPL are discretized within a radius R box . For this test we set R box = 20 fm. Since Supplementary Table 1.
The neutron emission widths Γ n and density of neutron vapor n gas for 238 U and 258 U as calculated with our model (FT-DIRHBz) using the relativistic DD-ME2 interaction, compared to results in Ref. [20].
Ref. [20] FT-DIRHBz interaction. In Fig. 8 we show the isotopic dependence of neutron emission lifetimes for samarium isotopes at T = 2 and 3 MeV. We observe that the agreement between the two calculation methods is excellent, all up to the neutron drip line. Furthermore, the influence of the particle continuum is especially pronounced at T = 3 MeV, reflected by more than an order of magnitude shorter neutron emission lifetimes for considered nuclei compared to T = 2 MeV. However, this example demonstrates that when supplemented with a proper continuum subtraction procedure, one can use the HO code to study neutron lifetimes within the scope of validity of Eq. (S29).
To further exemplify the validity of our approach, it is instructive to compare the results between different theoretical models. Here we compare our results with those from Ref. [20], which assume axially-deformed nuclei calculated with coordinate-space discretization and non-relativistic Skyrme SLy4 interaction. To this aim, we employ the axially-deformed FT-DIRHBz code with the DD-ME2 interaction and perform calculations for 238 U and 258 U at T = 1.0, 1.5 and 2 MeV. Our results are compared with corresponding results from Ref. [20] in Table 1. We observe that the neutron emission widths Γ n are consistent between both calculations, being of the same order of magnitude. For 238 U, as the temperature increases, the number of states in the continuum grows, thus increasing the neutron vapor density n gas . Following Eq. (S29), the density of neutron vapor states is proportional to neutron emission widths, and therefore, as the temperature increases, neutron emission lifetimes get shorter. Additional 20 neutrons in 258 U allow for more loosely-bound nucleons contributing to the neutron vapor thus increasing the neutron widths and decreasing the lifetimes for a fixed temperature.
We note the peculiarity in calculating the box radius R box within the axially-deformed HO code. The oscillator lengths are now represented by corresponding oscillator length in radial b ρ and z-direction b z , determined from the volume conservation relation b 3 Therefore, within the FT-DIRHBz, the discretization volume is still 4π/3R 3 box .

Supplementary Note 4. TEMPERATURE EVOLUTION OF QUADRUPOLE DE-FORMATION
Within the BLV subtraction procedure, the mean value of an observable ⟨O[ρ]⟩ T at temperature T , is a function of the subtracted density (ρ), defined as the difference between the density of the Nucl+Vap system (ρ) and Vap system (ρ). For the relativistic EDFs, the baryon density is equal to the vector densityρ v , which satisfies Eq. (S15). In the following, we present the results for the temperature evolution of isoscalar quadrupole deformation for selected even-even nuclei.
Starting from the proton(neutron) subtracted vector densityρ p(n) v the proton(neutron) quadrupole moment is defined as where (r ⊥ , z) are the cylindrical coordinates. It is more customary to express the results in terms of dimensionless variable β p(n) 2 defined as [21] β p(n) 2 where Z(N ) denotes the proton(neutron) number and R 0 = 1.2A 1/3 fm. The isoscalar quadrupole deformation is defined as β IS 2 = β p 2 + β n 2 . To demonstrate the effect of the shape phase-transition in more detail, we show the temperature evolution of β IS 2 for particular isotopes of neodymium (Z = 60) in Fig. 9(a)-(d), as calculated with three relativistic EDFs. For 150 Nd in Fig. 9(a) and DD-ME2 functional, we observe that at low temperatures, deformation can be slightly increased with increasing temperature. For T > 0.5 MeV, β IS 2 decreases smoothly with temperature up to T p c = 2 MeV, where it suddenly vanishes. The trends predicted with point-coupling functionals are similar, although T p c = 2.2 MeV for DD-PC1 and T p c = 2.4 MeV for DD-PCX. A higher phase-transition temperature is predicted with all functionals for 170 Nd in Fig. 9(b) where p c is close to 3 MeV. 180 Nd also has a prolate shape, which reduces to spherical at p c = 2.4 MeV for DD-ME2 , p c = 2.5 MeV for DD-PC1 and p c = 2.6 MeV for DD-PCX. For 190 Nd in Fig. 9(d), the shape phasetransition temperature is lower, being around 1 MeV for all functionals. However, the DD-ME2 predicts a slightly prolate shape of the nucleus, while both DD-PC1 and DD-PCX predict oblate configurations. Note that in the case of 190 Nd the ground-state difference between the oblate and prolate shapes is around 3 MeV, while for 150 Nd it is only around 0.1 MeV, signifying a strong configuration mixing.

Supplementary Note 5. SYSTEMATIC CALCULATIONS WITH DIFFERENT EDFS
To estimate the systematic uncertainties in predicting the two-nucleon drip lines at finite temperatures, it is instructive to compare results obtained using various EDF parametrizations. In Figure 10  interactions [22][23][24]. For instance, the meson-exchange DD-ME2 interaction consists of three meson-nucleon vertices in its mean-field part, corresponding to the exchange of the σ, ω, and ρ mesons, while the point-coupling interactions (DD-PC1 and DD-PCX) reduce the meson propagators to simple delta functions. Following Ref. [21], the systematic uncertainty is given as a spread across different EDFs. Assuming that A(X) is an observable of interest, being a function of different variables X, we define the systematic uncertainty of observable A as where A max(min) (X) represents the maximum(minimum) values of A(X) with different EDFs. Note that variables X here represent input parameters of the model, e.g., proton and neutron number or temperature T , X = (Z, N, T, . . .). We take the mean value as the average across different model predictions where the sum goes over the DD-ME2, DD-PC1, and DD-PCX interactions.
Statistical uncertainties stem from errors in optimizing the EDF parameters. For their calculation, one requires the so-called covariance matrix, in addition to the first derivative of the observable with respect to the model parameters [22]. In this work, we calculate the statistical errors for the recently optimized functional DD-PCX [9]. Assuming a linear behavior of the observable with small changes in model parameters in the vicinity of the χ 2 -optimum, we can calculate the first derivatives either by (i) linear regression or (ii) finiteelement derivatives. To this aim, we use the mesh as defined by the second-order Richardson extrapolation [23]. For a given step size h and observable A(p 0 ), where p 0 is the optimal parameter set, it requires evaluation of the observable at 6 points in the parameter space.